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ABSTRACT 







The scalar form of the approximate factorization method was used to develop a new code for the solution 
of three-dimensional internal laminar and turbulent compressible flows. The Navier-Stokes equations in 
their Reynolds-averaged form were iterated in time until a steady solution was reached. Evidence was given 
to the implicit and explicit artificial damping schemes that proved to be particularly efficient in speeding 
up convergence and enhancing the algorithm robustness. A conservative treatment of these terms at the 
domain boundaries was proposed in order to avoid undesired mass and/or momentum artificial fluxes. 
Turbulence effects were accounted for by the zero-equation Baldwin-Lomax turbulence model and the q-w 
two-equation model. For the first, an investigation on the model behavior in case of multiple boundaries 
was performed. The flow in a developing S-duct was then solved in laminar regime at a Reynolds number 
(Re) of 790 and in a turbulent regime at Re=40,000 by using the Baldwin-Lomax model. The Stanitz 
elbow was then solved by using an inviscid version of the same code at M jn | et =0.4. Grid dependence and 
convergence rate were investigated showing that for this solver the implicit damping scheme may play a 
critical role for convergence characteristics. The same flow at Re=2.5-10 6 was solved with the Baldwin- 
Lomax and the q-u> models. Both approaches showed satisfactory agreement with experiments, although 
the q-w model was slightly more accurate. 


INTRODUCTION 

With the advent of more and more powerful supercomputers, the numerical solution of three- 
dimensional turbulent flows became possible (ref. 1). Although it is well known that lower order turbulence 
closures fail to reproduce secondary motions of Prandtl’s second kind (turbulence driven), they normally 
succeed in predicting secondary flows of the first kind (pressure driven). Therefore, for many complex 
configurations it has been possible to obtain fairly accurate results with zero- and two-equation turbulence 
models (ref. 2) with a reasonable prediction of pressure driven secondary flows. 

For three-dimensional blade-passage flows, a correct prediction of the wake behavior was obtained by 
Yokota (ref. 3) by means of the standard high-Reynolds-number form of the k-e two-equation model with 
the wall function approach. However, the zero-equation model implemented in reference 3 did not give a 
satisfactory depiction of the flow in the wake region. For the prediction of blade pressure distribution, the 
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Baldwin-Lomax (ref. 4) zero-equation turbulence model was found to give accurate results in several flow 
configurations (ref. 5) despite the low convergence rate. For incompressible internal flows with no 
separation, good results were obtained by Towne (ref. 6) with a zero-equation turbulence model and a 
parabolized Navier-Stokes solver. Still, for practical flow configurations, quite long computing times are 
necessary mainly because of the large number of points usually required for a detailed description of the 
flow field. Moreover, the nonlinearities associated with nearly all turbulence models can play a significant 
role in slowing down the convergence rate to the steady state solution. This behavior occurs in both the 
implicit and explicit flow solvers and is intrinsic to the turbulence models, zero or two-equation. 

In two dimensions, a wide variety of flow conditions have been accurately solved by means of low- 
Reynolds-number forms of the k-e model in which the effect of laminar viscosity is explicitly accounted for 
(refs. 7 to 9). In nearly all of the flow conditions investigated, these forms proved to be more accurate than 
the standard high-Re formulation, provided that a sufficient number of grid points were located inside the 
viscous and buffer layers; in fact, secondary motions and losses are mainly driven by what happens close to 
walls so that a correct description of this flow region is crucial for an accurate simulation of the flow 
pattern. Rodi (ref. 9) found that the low Reynolds number forms of the k-e model could predict secondary 
flows that are normally lost with the high Reynolds number form. Unfortunately, the first author (ref. 7) 
found that some of these forms were extremely stiff from a numerical point of view. The stiffness was 
mainly caused by the low-Reynolds-number effect terms in which exponential functions are introduced to 
model the wall effects. 

From this standpoint, it appeared worthwhile to investigate some features of turbulence models for 
internal turbulent flows by using an implicit algorithm. Since complex flow patterns, such as separation 
and dominant viscous effects, are expected in internal flows, the implicit approach was selected to increase 
the robustness and convergence rate of the numerical procedure when zero and two-equation turbulence 
models are used. 


DESCRIPTION OF THE ALGORITHM 

Navier-Stokes Equations 

The Boussinesq hypothesis allows the turbulent shear stresses to be related to the mean strains via the so- 
called eddy viscosity so that, under this assumption, the three-dimensional Reynolds-averaged compressible 
Navier-Stokes (N-S) equations can be written in divergence form and, subsequently, be transformed from 
the Cartesian coordinate system (x,y,z) to the generalized curvilinear coordinate system (£,??, C)- The 
resulting set of equations can be written in vector form as follows: 

dQ , <9E , 8F , dG _ <9E V , dF v , <9G V 

dt + d(, + dr] + 8( ~ d£ + drj + d( 

where Q is the vector of unknowns; E, F, and G are the flux vectors of the convective terms; E v , F v , and 
G v are the flux vectors of the diffusive terms, the complete definition of which can be found in reference 19. 

In the set of equations, p is the fluid density, p is the static pressure, e is the total energy per unit 
volume, a is the sound speed, and (U,V,W) are the Cartesian components of the velocity vector. According 
to the Boussinesq assumption, the diffusion coefficients for momentum and energy are defined as follows: 

/"eft = A*l + Ht and 

U t 


in which /q is the laminar viscosity, which was considered to be independent of the static temperature, and 
p t is the turbulent viscosity obtained from the turbulence model. In this set of calculations the turbulent 
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Prandtl number Pr t was set equal to 0.90 and the laminar Prandtl number Pr| was 0.72. 

Discretization 

The flux vectors are discretized by using centered finite differences. The metrics are usually obtained from a 
chain rule expansion of x^, x^, x^, y>, y^, y^, z^, Zj j, and z >. When the centered discretization is used for 
the metrics, it can be shown that in three dimensions the metric invariants are not satisfied (ref. 10). This 
may result in a large discretization error. However, it is possible to satisfy the invariants by a simple 
averaging technique that gives metrics that are similar to those computed by a finite volume method. For 
example, £ x is computed as follows 

in which J is the Jacobian of the coordinate transformation and \ is a central average operator. This 
averaging process was found to ensure better mass conservation properties in the present calculations 
expecially for highly stretched grids. 

Approximate Factorization; Scalar Form 

The approximate factorization method first proposed by Beam and Warming (ref. 11) splits an n- 
dimensional operator into the product of n one-dimensional operators. This technique provides a strong link 
between the equations insofar as they are solved fully coupled. The main drawback of this method lies in 
the necessity of a time-consuming block tridiagonal or pentadiagonal matrix inversion. This problem 
becomes more evident in three dimensions where the coupled solution yields a 5x5 block tridiagonal 
matrix. In order to make this algorithm more efficient and still maintain its strong implicit nature Pulliam 
(refs. 10 and 12) proposed a scalar form of the approximate factorization. The form of the standard 
algorithm in three dimensions can be written as follows: 

[ I+0At(<yA - <^A V )]*[ I+0At(fyB - <5^B V )]*[ I+GAt^C - $JC V )]*AQ=RHS 

in which I is the identity matrix; 0 is a parameter that allows the explicit-implicit nature of the space 
operator to be weighted (in the present calculations, since the steady state solution was sought, 0=1 was 
used which gave first-order accuracy in time); 6 is a centered difference operator; At is the time step; 
AQ=Q n+1 - Q n ; Q=J Q; A, B, and C are the convective Jacobians; A v , B v , and C v are the diffusive 
Jacobians in the three directions £, 77 , and (; and RHS represents the convective and diffusive fluxes and is 
defined as 


RHS=At( 


QM.dF.dG 

d£ dr) d( 


<9E V 

d* 


, dF v , <9G V 
+ ~d^ + ~dc 


) 


The Jacobian matrices A, B, C, A v , B v , and C v , which can be found in Pulliam (ref. 10), are full and 
require a general inversion process. Pulliam and Chaussee (ref. 12) proposed a more efficient procedure by 
making some simplifications. First, the hyperbolic property of the convective Jacobians allows the 
diagonalization as 


A = T, A^ Ty 1 ; B=T 71 A tj T y 1 ; C = T ( A ( Ty 1 

in which T are the eigenvector matrices, defined in reference 12. The eigenvalues of the Jacobians in three 
dimensions, A, are given by 
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( 3 ) 


A £ = D[U,U, U, f/-aJ^ x 2 +^ y 2 +G 2 ] 

A V = D[V,V,V, F+a \ V x. 2 +Vy 2 +Vz 2 , F-aJ r ]x 2 +r 1y 2 + Vz 2 } 

A < = D[W,W, W, VF+aJ Cx 2 +Cy 2 +C Z 2 , W- a J C X 2 +Cy 2 +Cz 2 ] 

where U, V , and IF are the unsealed contravariant velocities and D stands for main diagonal only. It is 
now possible to introduce equations (2) into equation (1) with eigenvalues defined in equation (3). Doing 
so, it is impossible to diagonalize both the convective and diffusive Jacobians since they have a totally 
different set of eigenvectors. Neglecting the diffusive Jacobians and assuming the eigenmatrix to be locally 
constant, equation (1) can be rewritten as 


T^[I+0At^A^]*N*[H-0At^ 7? A^j]*P*[I+0At^A^]=KT^ 1 *AQ=RHS 

( 4 ) 

in which the two matrices, N=T^ 1 and P=T^ T> have the nice property of being solution independent 
so that they can be computed only once (ref. 12). Obviously, equation (4) is only an approximation of the 
full form but it requires only a scalar tridiagonal or pentadiagonal matrix inversion since the A matrices 
are diagonal; this implies a reduction of nearly 50% of the operations required by the standard algorithm. 
We note that, even though the interior domain is solved implicitly, the boundary conditions are imposed in 
a fully explicit manner. This was done by setting AQ=0 at the domain boundaries during the implicit 
sweeps. 

Since the steady state solution is sought, the following local time-stepping formula, based on the 
approximate constant CFL condition is introduced, in which the contribution of the diffusive terms is 
accounted for 


At^ — \U\ + a .j^x + 6y + 6z + f'eff 1 (£x + £ y + 6z) 

At,; = W\ + a ^ + rfi + 77 f + p eff Re' 1 (77I + rfi + rfc) 

At^ = \W\ + a JCx + Cy + Cz + /^eff ^- e 1 (Cx + Cy + Cz) 

At = CFL 

At ^ -{- At rq At q 


Implicit Treatment of Diffusive Terms 

For inviscid solutions, the only assumption done to derive the algorithm is that the differentiation of the 
eigenmatrices is neglected in equation (4). In addition, when solving the N-S equations, the Jacobians of 
the diffusive terms are normally assumed to be negligible. This does not cause any stability or convergence 
problems for external flows where the diffusion dominated region is small and restricted to a limited part of 
the computational domain. But this simplification may cause troubles for internal flows where the diffusion 
dominated region can be large. According to this, for internal flows it is convenient to introduce an 
approximation of the eigenvalues of the diffusive terms jacobian and put it into equation (4). Two forms of 
this approximation have been considered. 

0 Pulliam’s approximation. 

Pulliam (ref. 10) proposes an approximate form of the diffusive Jacobian eigenvalues that was obtained by 
intense numerical testing. This form is 
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( 5 ) 


A* = {p p eff Re ' 1 (£ x 2 +£y 2 +£ z 2 ) r 1 ) ■ D [1,1, 1,1,1] 

= (p Ai eff Re ' 1 (% 2 +7? y 2 +^ 2 ) J" 1 ) • D [1,1,1, 1,1] 

A* = (p p eff Re " 1 (Cx 2 +Cy 2 +Cz 2 ) J" 1 ) • D [1,1, 1,1,1] 

In the present set of calculations it was found convenient not to weight the eigenvalues with the Jacobian 
of the coordinate transformation J -1 . These expressions are included on the implicit side of equation (4): for 
instance, the £ direction implicit operator reads 


[I + 0 At (6^ A^ - ] 


• Present approximation. 

The exact form of the diffusive terms Jacobian can be computed from the related flux vectors. For the 
three sweeps the main diagonal of such matrix may be conveniently approximated as 


A £ V = D [0,£*£, «£,<*£, 7 Pr" 1 ^] 

A ij = D [0,a r j,a r j,a 1 rj,j Pr 1 a^] ( 6 ) 

A^ v = D [0 , 7 Pr" 1 ^] 


in which 

a £ = Peff R e_1 J " 1 (£x 2 +£y 2 +£z 2 ) 
a 7j = Meff Re' 1 J' 1 (Vx 2 +Vy 2 +Vz 2 ) 
«£ = /i eff Re -1 J ’ 1 (Cx 2 +Cy 2 +Cz 2 ) 


Regarding the extra-diagonal terms as negligible, the previous diagonal matrices are a good approximation 
of the diffusive terms Jacobians and can be put into the implicit side of equation (4): for instance, the £ 
direction implicit operator reads 


[ I + 0 At (S ( A f - ^ Ap ] 


in which the diffusive terms contribution is approximated as a first order derivative. It is possible to prove 
that the first approximation increases the main diagonal dominance by summing to it an extra term, 
whereas in the second approximation an artificial term is subtracted from the off-diagonal components 
while leaving the main diagonal unchanged. A comparison of the two approaches, equations (5) and ( 6 ), 
was performed on a simple straight channel geometry at Re~1000 and M jn | et =0.3 in a laminar flow 
regime. For this very simple flow configuration there were no differences in convergence rate between the 
two approaches, and it was also possible to drop the diffusive terms on the implicit side of equation ( 4 ) 
without altering convergence. Differences started to appear at Re=50 because of the highly diffusive nature 
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of the flow. Figure 1(a) shows the best convergence history of the algorithm without any implicit treatment 
of the diffusive terms that were obtained at CFL=5. (The lower curve refers to the averaged residuals, and 
the upper curve refers to the maximum.) Figures 1(b) and 1(c) refer to equations (5) and (6). There are no 
appreciable differences between the two convergence rates obtained at CFL=10 since both curves show 
nearly the same slope. The same result could be obtained if the Baldwin-Lomax turbulence model was 
used. At least for this class of flows the approximate implicit treatment of diffusive terms given by 
equation (6) did not prove to be more efficient than equation (5). Further testing is necessary in order to 
verify this result at higher Mach numbers when differences in the convergence between equations (5) and 
(6) may appear since the density derivatives do not tend to vanish for compressible flows rather than for 
incompressible flows. 

(a) - no. diff. terms (i) - approximation (5) (c) - approximation (6) 



ITERATION N. ITERATION N. ITERATION N. 

Figure 1. Convergence tests for implicit treatment of diffusive terms. 


Added Explicit and Implicit Damping 

When using centered finite differences, it is necessary to introduce an artificial damping scheme in order to 
prevent odd-even velocity-pressure decoupling that occurs whenever the local Peclet number exceeds 2 in 
large-gradient regions. Moreover, the three-dimensional approximate factorization method can be proved to 
be unconditionally unstable since, when performing a linear stability analysis, one of the amplification 
factors is found to be slightly larger than one. This weak instability can be easily overcome by using an 
artificial damping scheme. From the two-dimensional version of the approximate factorization, it is also 
known that fourth order artificial damping is necessary to damp the numerical modes associated with the 
highest frequencies. Following the basic guideline given by Jameson et al. (ref. 13) and Pulliam (ref. 14), a 
blended implicit-explicit second plus fourth order nonlinear damping scheme was introduced in the present 
solver. 

© Explicit 2 nd and 4 th order damping. 

In the original formulation, the artificial damping is simply equally scaled according to the local At in the 
three space directions. This could be done mainly because the damping scheme has been first applied to 
in viscid flows. For viscous flow calculations, Pulliam (ref. 10) found it convenient to scale the damping 
terms according to the directional spectral radius u. The scheme for the £ direction sweep yields 
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( 7 ) 


d^ 2) +d^ 4) 


^(((aJ- 1 ) i+1Jik +(^J- 1 ) iJik )(a ) fJ ) k ^^ (Q |Jik ) - wfj.k^QiJ.lt)) 


(2) (4) 

The complete definitions of u and w are given following Pulliam 


T — 
i , j , k 


Pj+1 j,k 2 Pj,j,k Pj-l,j,k 
Pi+lJ.k + 2 Pi,j,k + Pi-1, j,k 


“u’k = £i<2) it .«<T UJi( , T |Jik , T i+1Jik ) 


U) 


(4) 


j k = max( 0., At - wj'j k ) 


.( 2 ) 


( 8 ) 


and the unsealed spectral radius is defined as 

^i.j.k = 1^1 + a Jd + £y + £z 

Since the present set of calculations has been performed for shock-free flow fields, the shock sensor 
term defined in equation (8) was switched off by always taking Yj j k . Nevertheless, this modification does 
not really affect the damping scheme since the second order term is active only in presence of shocks. 
Actually, this scheme gives only fourth order damping in smooth shock-free flow fields so that in most of 
the calculations fr was set equal to zero whereas the fourth order weight, fr , ranges from 1/16 to 
1 / 12 . 8 . 

Equation (7) is then added to the RHS of equation (4), 

RHS = RHS + At (D^ 2) + D^° + D^ 2) + D^ 4) + D^ 2) + D^ 4) ) 


® Implicit 2 nd or 4 th order damping. 

To enhance the algorithm stability and convergence rate, it is helpful to include the artificial damping 
terms on the implicit side. For both the 2 nd and the 4 th order damping, the implicit treatment augments 
the diagonal dominance of the scalar system with a beneficial effect on the convergence rate. As it was 
pointed out by Pulliam (ref. 10), there is a stability limit for the weights that can be used for the artificial 
damping terms. This limit is actually connected to the magnitude of the amplification factor of the scheme 
modified with the artificial terms. To obtain the best convergence rates the implicit damping had to be the 
exact Jacobian of the explicit counterpart added to the right-hand-side. In case the second order damping is 
treated implicitly, equation (4) must be modified to include the added implicit terms: for instance, the £ 
sweep implicit side will be 


{‘+6At(« £ A { - « £ (((rf- 1 ) i+1 j ik +(rf- 1 ) lijik )(u,fJ > R 6 £ J iJ , k 



(») 

This form maintains the tridiagonal nature of the Jacobian matrix. 

When treating the fourth order damping implicitly, the differentiation of the Jacobian matrix brings 
about a scalar pentadiagonal system that can be written for the £ sweep as 


{l+0At(5^ - ^(((^J- 1 ) i+1Jik +(^J- 1 ) kjik )( W [ 4) k ^J iijik )))} 


( 10 ) 
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A brief set of tests on a straight channel geometry proved that the 4 th order implicit damping (eq. (10)) 
could bring a large gain in convergence rate with respect to the 2 nd order (eq. (9)). For a straight three- 
dimensional channel with approximately 10 4 points with an inlet-section-width to channel-length ratio of 
15 (typical of internal flows geometries), the tridiagonal solver associated with equation (9) could be run at 
CFL=2 which gave the convergence history shown in figure 2(a), whereas the pentadiagonal solver 
associated with equation (10) proved to be much more robust and gave the much higher convergence rate 
shown in figure 2(b) with CFL=10. Despite an increase of around 25% in computational time to invert the 
pentadiagonal matrix with respect to the tridiagonal one, the implicit 4 th order option proved to be much 
faster and more robust, and was retained in all of the calculations. 


(a) tridiagonal 


( b ) pentadiagonal 



Figure 2. Convergence tests for tridiagonal or pentadiagonal solvers. 


@ Wall treatment of artificial damping. 

It is well known that any kind of artificial damping term introduces an error that has to be minimized in 
order not to affect the solution heavily. Among the various approaches that can be found in the literature, 
the nonlinear damping formulation proposed by Jameson et al. (ref. 13) ensures that the second order 
terms are introduced only at shocks while keeping the fourth order in the smooth region of the 
computational domain. However, it is possible to prove that, if the presence of the boundaries is not 
properly accounted for when introducing the artificial terms, nonzero momentum and mass fluxes can be 
produced at the boundaries. This fact can be easily seen by considering figure 3. Figure 3(a) shows the 
fourth order damping weights applied at every single point i for the £ direction. The fourth order difference 
stencil used here is 


^ ( Q i j , k ) — Qi-2,j,k~4 Qj-l,j,k+ 6 Qi,j,k‘4 Qi+I,j,k + Qi+2,j,k 


The fourth order damping is normally switched off at i= 1 and i= 2 so that no special treatment at the wall 
is required. Figure 3(a) shows that, doing so, the sum of the fourth order damping weights for every 
locaton i is zero only in the internal flow field for i> 5. This ensures that no net fluxes are added only in 
the internal domain, while the sum of the artificial damping weights yields nonzero values for l<i<4. 
These weights actually correspond to a nonzero third order derivative centered at i=(2+l/2). The same 
procedure may be followed for the second order damping that is switched off at i= 1 (figure 3(b)). The sum 
of the weights for every location i is zero only for i> 3. This yields a first order derivative centered at 
i=l/2, that corresponds to a first order flux at the wall. These flux errors can be easily controlled in two 
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dimensions by a grid refinement in the wall proximity where high gradients are expected. In three- 
dimensional internal flows we are forced to use coarser grids and the wall boundary condition is applied on 
very large surfaces with the result of possible large mass errors. To control the mentioned flux errors, a 
third order derivative, with the differencing stencil given by 

^(QiJ.k) = -Qi-lj,k+ 3 Qi,j,k' 3 Qi+1 j,k+ Qi+2,j,k 

was added at i— 2 to balance the fourth order damping weight. The same procedure was followed to balance 
the second order damping at the wall where a first order derivative, with the differencing stencil given by 


was added at i= 2. 


^(Qi,j,k) — Qi-lj,k + Qi,j,k 
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b) SECOND ORDER 


Figure 3. Artificial dissipation treatment at boundaries. 


TURBULENCE MODELS 


Baldwin-Lomax for Multiple Boundaries 

The standard version of the Baldwin-Lomax (ref. 4) zero-equation turbulence model was implemented here. 
This two-layer model divides the flow field into an inner layer close to the wall, in which viscous effects are 
dominant, and an outer layer. The influence of the solid wall is damped according to the Van Driest 
exponential function. 

This algebraic turbulence model was originally developed for a single boundary layer, but in three- 
dimensional internal flows the presence of multiple boundaries causes the interaction of more than one 
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boundary layer. Whereas the inner, or viscous layer, is driven by what happens on the closest wall only, for 
the outer layer an averaging procedure is necessary to account for the various wall effects. In the present set 
of calculations three approaches have been examined to account for the presence of more than one 
boundary layer. 

e Wall Treatment 1 

Only the geometrically closest wall is considered to compute the outer viscosity without any averaging. 

© Wall Treatment 2 

The inner layers are driven by the closest wall only, whereas the viscosity in the outer layer is computed 
with a simple weighted average according to the inverse of each wall distance as follows 


^t, outer 


n ,, w 

1 /^t , outer 

W w j 

VWw 


where W w is the wall distance, w is the wall number, and N is number of walls present in a cross section. 

© Wall Treatment 3 

In case only one boundary layer is present, the Van Driest damping succeeds in modeling the wall effect. 
Starting from this standpoint, the inner layer viscosities are computed using only the closest wall 
contribution, while the viscosity in the outer layer is computed as a simple weighted average according to 
the inverse of the value of the Van Driest damping expression for each wall. 

q-w two-equation Model 

In a previous investigation, Michelassi (ref. 7) found that the low-Reynolds number forms of the two- 
equations models, like k-e, could give accurate prediction of two-dimensional incompressible separated 
flows. Unfortunately, these forms were found to be numerically stiff, mainly on account of the correction 
terms introduced to model the low-Re effects that necessitate a strong mesh refinement in the sublayer. 
Furthermore, an initial profile for the turbulent quantities has to be specified consistently to start the 
calculations. A first attempt to implement the Chien’s and the Rodi’s two-layer low-Reynolds number 
forms of the k-e model did not bring any converged result mainly because of difficulties in specifying both a 
sufficient mesh refinement and proper initial profiles for complex three-dimensional flows. Coakley (ref. 15), 
reassembling the Jones and Launder low-Reynolds-number form of k-e model, proposed the q-o> two- 
equation model, in which the effect of molecular viscosity is directly modelled: This formulation ensured 
better numerical behavior as compared with other low-Re formulations. This vector form of the model 
transport equations rewritten in our notation is (ref. 19) 


dQ 

dt 



<9F | 8G dEv , dF v 

dr] d( d£ drj 


, <9C V , n 
+ ^- + H 


( 11 ) 


in which H stands for the sink and source terms vector. The full definition of the flux vectors can be found 
in reference 19. The two transported quantities, q and u, are related to the more familiar k and e via the 
following relations: 

q = k 1 ^ 2 ; w = e / k (12) 


It is important to observe that e appearing in equation (12) is the isotropic part of the dissipation 
rate; this quantity does not account for any nonisotropic effect (for example, the presence of a wall) and 
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tends to zero on solid boundaries. (Conversely, the total dissipation rate tends to a finite value related to 
the wall shear stress.) This choice allows w to be used as an unknown since, assuming that both k and e are 
going to zero at the wall with the same slope, w tends to a finite value. 


q-w Solution 

The two transport equations for q and w are implicitly solved by using the same algorithm given in 
equation (1). The two equations are solved in a sequential manner and decoupled from the flow variables 
mainly because the coupling is provided only by the coefficients of the diffusive terms and the sink and 
source terms. Because of this choice, the scalar three-diagonal algorithm was implemented for the 
turbulence model solution. The only difference with respect to the solution of the N-S equations is the 
presence of the sink-source vector H in equation (11). This term can be included in the three sweeps: 


[ I+0At(-Hj c^+ S^A - <5|A V )]* 

[ I+0At(-Hj Cr) + S^B - ^B v )]* 

[ I+0At(-Hj C( . + S^G - <5);C V )]*AQ = RHS 


(13) 


where the same definitions given for equation (1) hold with the only addition to the Jacobian being the 
sink and source terms, Hj, that are weighted in the three sweeps according to c^, c^, c^, and c^+c^+c^.= 
1. This Jacobian is computed neglecting the contribution of the damping function D. Its form for the two q 
and w equations is 


H; 


d H q 

5(pr x q) 


1 Cp DS 

0 U) 


p c o 

id ~ 7T 


Hf 


U> 


<9H 
d(pT 1 w) 


1 Ci D P d 


Co 2 w 


in which S and P d are production terms (ref. 19). In place of their exact form, Coakley (ref. 15) proposes 
an approximation of the Jacobians based on the turbulent viscosities that should ensure the dominance of 
the main diagonal. Figure 4 shows the comparison of the convergence rates of every single variable 
obtained without any sink or source terms Jacobian, with the exact Jacobian, and with Coakley’s 
approximation in a typical internal flow geometry. Surprisingly, there is no big gain in introducing the 
Jacobian in the implicit side of the operator since with the three formulations it was always possible to 
obtain the same residual reduction in 300 iterations. The choice of the sweeps in which the two Jacobians 
H? and Hj*' are introduced is not important. The error introduced in the approximate factorization of the 
implicit side of equation (13) increases roughly a factor proportional to Hj when introducing the Jacobian 
in the three sweeps, thereby choosing c^_=c T j=c^= 1/3: this has only a weak influence on the convergence 
rate. Nevertheless, in the present calculations typically the best convergence rates have been obtained by 
using c.£=0, Crj=Q.5, and 0 ^= 0 . 5, where £ is the main flow direction and rj and £ are the fine grid 
directions. 

Although physical evidence shows that the turbulent kinetic energy k is zero at solid walls, the 
boundary condition for e, and consequently w, is less evident. For the q-w model, Coakley (ref. 15) found it 
convenient to impose a zero-normal derivative at the wall; this condition was retained in the present 
calculations. 
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Figure 4. Convergence tests for implicit treatment of q-w sink-source terms. 
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Incompressible S-duct 

A first validation of the code was performed by computing the flow in an S-duct with a constant area cross 
section. The measurements (ref. 16) were taken for an incompressible fluid (water). The S-duct is given by 
two 22.5° bends with 4-cm hydraulic diameter and 28-cm mean radius of curvature. This geometry was 
regarded as an interesting test since flow passages with similar shapes are often used to redirect the flow for 
air intakes in aeronautical engines. The efficiency of such ducts may be heavily affected by the presence of 
secondary flows so that the ability of a code to detect secondary velocities is essential for a proper depiction 
of the flow pattern. 

The sketch of the flow domain is shown in figure 5. Since the geometry and the flow were symmetric 
with respect to the x-r plane, it was possible to study only one-half of the duct by imposing a symmetric 
boundary condition on the same plane. For the laminar flow regime, the Reynolds number, Re, based on 
the inlet bulk velocity, Ub, and hydraulic radius was 790. The grids employed for this calculation are 
shown in figure 6. The coarser grid (figure 6(a)) has 70x29x15 points with a ratio between two 
consecutive grid cells in the cross flow directions of 1.1, and the refined one (figure 6(b)) has 80x59x30 
points with the same stretching ratio in the cross flow plane. The fluid adopted in the experiments is water 
so that, in order to have negligible compressibility effects, the isentropic Mach number was set equal to 0.1 
(which gives the inlet total pressure, outlet static pressure ratio) to ensure minor density changes. 

The flow pattern and the growth of the secondary motions is mainly pressure driven because of the 
very smoothly bending walls that cause no flow separation. Still, this kind of flow necessitates a very 
accurate prediction of the boundary layer, even in laminar regime, otherwise the secondary flows may be 
incorrectly predicted or completely lost. Figure 7(a) shows the measured (circles) and computed velocity 
profiles with the two grids mentioned above (dashed line is the coarser grid; solid line is the refined grid) 
for the 5 sections (indicated in fig. 5) on the symmetry plane reported in the experiments. The slight 
discrepancies in the computed velocities with the two grids at section 1, where the double-S starts, may be 
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attributed to the thinner boundary layer predicted with the refined grid that produced a flatter velocity 
profile. Nevertheless, the agreement with experiments is fairly good. The agreement does not deteriorate 
for sections 2, 3, and 4 in both the coarse and refined grids. Section 5 clearly shows that the secondary 
motions could not be predicted by the coarse grid, whereas they are fairly well reproduced by the refined 
one. This behavior is thought to be independent of the grid points in the main flow direction, (as it will be 
demonstrated in a further test) in which only 10 points are added in the refined grid, and it is closely 
related to the poor resolution of the boundary layer provided by the coarse grid in which the cross-flow 
momentum diffusion is clearly overestimated. 

Figure 7(b) shows the velocity profiles for the midspan section (r=l/2) at the same five sections. 
Basically the same comments given for the symmetry plane could be repeated here for the first three 
sections, whereas on section 4 and 5 agreement deteriorates in proximity of the symmetry plane. This may 
be attributed to the poor grid quality close to the symmetry plane, but also to the fact that in the 
calculations a zero gradient condition was imposed in the direction normal to that plane, whereas 
experiments show that the velocity gradient is far from being zero on the same plane for section 4. The 
typical computed secondary motion pattern at the exit of the second bend is shown in figure 7(c): The 
complex flow pattern exhibits two counter-rotating vortices in proximity of the two curved walls. 

Nearly 1500 iterations were necessary to reach an averaged residual of the order of 10 5 on the 
refined grid, but typically less than half of these iterations were necessary to obtain the same residual with 
the coarser grid. These calculations were performed with an early version of the code where only the second 
order implicit damping was implemented so that the slow convergence rate may be attributed to both the 
small Mach number and to the small CFL number that could not exceed 2 without encountering stability 
problems. 

The turbulent flow regime was run at the experimental condition of Re=40.000. For this calculation, 
it was necessary to provide a more stretched grid at the wall, so that the 80x59x30 grid was reassembled 
with a point expansion ratio in the cross flow directions equal to 1.3. This provided the necessary point 
clustering at the boundaries to describe the thin boundary layer and allowed the first point at the wall to 
be placed at y’*’w2. For this test case, the Baldwin-Lomax model was implemented with wall treatment 1 
(see section on Baldwin-Lomax for multiple boundaries). 

Figure 8(a) shows the set of measured and computed velocity profiles on the symmetry plane for the 
same five cross sections given in figure 5. For this flow configuration the agreement is reasonably good, 
especially at cross sections 4 and 5 where the pressure gradients induce a strong secondary motion that 
seems to be correctly predicted by the present solver insofar as the agreement with measurements does not 
deteriorate as the second bend exit is approached and the secondary velocities reach their maximum. The 
momentum transfer to the external part of the second bend, particularly evident in sections 4 and 5, is 
reproduced well since the computed velocity profile asymmetry seems to be in close agreement with 
experiments. 

Figure 8(b) shows the computed and measured velocities on the midspan plane along the duct. The 
agreement is again good for the five sections and better than that found for the laminar flow regime mainly 
because the measured velocity profiles exhibit a zero gradient on the symmetry plane that is correctly 
modelled by the symmetry condition imposed on the x-r plane. Still, at sections 3 and 4 the kink in the 
velocity profile close to the wall is not correctly predicted. 

The static pressure coefficient in the main flow direction, defined as 


C p = 


P ' Pinlet 



is shown figure in 8(c) at three different locations (z=0,r=0), (z=0,r=l), and (r=l/2,z=l/2). The 
agreement is generally good. The pressure trend is correctly predicted together with the head loss. 
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Figure 8(d) shows the convergence history obtained with CFL=2: For this calculation the same early 
version of the code mentioned previously was used. The solid upper line refers to the maximum residual, 
and the dashed one refers to the averaged residual; the spikes present in the first curve correspond to the 
updating of the turbulent viscosity performed every five iterations. It is remarkable that the pressure 
distribution did not change after the first 500 iterations, whereas to get the correct velocity profiles, it was 
found necessary to reach a residual of the order of 10" 5 to 10~ 6 . 


Stanitz Elbow 

The flow in the accelerating rectangular elbow with 90“ turning and variable cross section described by 
Stanitz (ref. 17) has been computed with an inviscid version of the code and with both the Baldwin-Lomax 
and q-w turbulence models. This test case was selected since it provides a good set of measurements 
including wall pressure distribution and visualization of secondary flows at the elbow exit section. The 
shape of the elbow was analytically computed by Stanitz to give no separation with a strong area reduction 
and a specified pressure distribution on the side wall under incompressible flow conditions. A sketch of the 
experimental setup is shown in figure 9. The flow and the elbow geometry are symmetric with respect to 
the x-y midspan plane thereby allowing a zero normal gradient condition. Among the various flow 
conditions investigated in reference 17, we selected the one with M exjt =0.4 and with no spoiler at the duct 
inlet with a thin initial boundary layer. 

© Inviscid Calculations 

A first set of tests was performed with an inviscid version of the code. The convergence characteristics 
of the scalar form of the approximate factorization could be tested for inviscid calculations where the 
necessity of accounting for the diffusive terms on the implicit side of equation (4) drops. Only the pressure 
distribution on the walls of the elbow was compared with measurements in this set of calculations: no 
attempt was made to specify an experimental inlet profile of total pressure that was kept flat. The inlet 
boundary condition was specified extrapolating the Riemann invariant from the first section inside the duct 
(ref. 18). Two grids were used: a 25x15x11 grid and a 50x15x11 grid with constant grid spacing in the 
cross-flow directions. The grid point locations in the streamwise direction were made to coincide with the 
points supplied by Stanitz (ref. 17) for the description of the elbow geometry with the addition of three 
cross sections at the outlet to allow the use of a zero gradient condition. 

Figure 10 shows the qualitative static pressure and Mach number isolines on the symmetry plane of 
the elbow. The small wiggles visible at the domain exit are due to the very small 4 th order damping weight 
that was set equal to fir =1/256, together with fir =0. The small value of the fourth order damping 
weight is allowed by the very coarse grid in the cross flow direction that automatically introduced a 
numerical diffusion. Despite the very coarse grids implemented here, the static pressure distribution p s , 
defined as 


P = P Pexit 
s Ptotal ' Pexit 


shown in figure 11, reveals a fairly good agreement with experiments. Although the pressure drop position 
is located correctly for the two grids on both the side and symmetry planes, the kink in the static pressure 
distribution on the suction surface at the side wall could not be reproduced since it is caused by the 
presence of secondary flows. The pressure distribution appears to be independent of secondary flows induced 
by viscous effects until the first part of the bend is reached. The local pressure rise located at 5=2 (where S 
is the streamwise coordinate along the centerline) is introduced by the growth of secondary velocities and is 
totally lost by the inviscid solver. The solution proved also to be fairly grid independent, at least to the 
grid refinement in the main flow direction: No additional tests were performed to verify the influence of the 
point distribution in the cross-stream direction. 
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The implicit treatment of 2 nd and 4 th order artificial damping was compared using the 50x15x11 
grid. Figure 12 shows the comparison between the two implicit damping schemes; the solid line refers to 
the 2 nd order implicit damping with Q =1/4, and the dashed line refers to the 4 th order implicit option 
with fr =1/256. The gain in convergence rate is remarkable; in fact, the fourth order solver could be run 
at CFL=10, whereas the second order one could not be run at CFL>5. Any further increase of the artificial 
damping weights gave slower convergence histories. Unfortunately the same convergence rate is not 
obtainable for viscous calculations because of the strong point clustering and viscous effects not exactly 
accounted for on the implicit side of the operator. 

© Viscous Calculations 

The viscous calculations in turbulent flow regime were performed on the five grids summarized in 
table I. The use of various point clustering and distribution allowed a comprehensive investigation on the 
mesh dependence of the calculations. With this set of grids it was possible to verify the influence of the grid 
points number in the main flow direction, with 51 or 99 points, and the cross flow direction with 31x21 
and 41x31 points with expansion ratios of 1.2 and 1.3. The refined grid (number 5) shown in figure 13 
adopts the same distribution of points in the main flow direction that was used for the refined grid in the 
inviscid calculations and which allowed placing the first grid point at the wall at y^«l. The Reynolds 
number, Re, based on the total conditions at the inlet section is approximately 2.5 - 10" 1- . 


grid 

number 

points 

expansion 

ratio 

1 

50x31x21 

1.2 

2 

50x41x31 

1.2 

3 

99x31x21 

1.2 

4 

50x31x21 

1.3 

5 

51x41x31 

1.3 


Table 1. Grids for viscous calculations on Stanitz elbow. 


These calculations were mainly aimed at the proper prediction of the wall pressure distribution that is 
heavily affected by the growth of secondary velocities and a qualitative comparison of the secondary flows 
predicted by the Baldwin-Lomax zero-equation model and the q -u> two equation model. The choice of the 
experimental spoilerless configuration allowed the turbulence models to be compared for a very thin 
boundary layer that required a heavy point stretching at the wall. Regarding the inlet boundary condition, 
in the Baldwin-Lomax model the inlet turbulent viscosity was extrapolated from inside the domain, and in 
the q-LO model a flat turbulent kinetic energy profile with various turbulence levels was specified at the inlet 
section while u> was extrapolated from inside the domain. 

The first set of tests concerns the Baldwin-Lomax model with the different wall treatments 
mentioned. With the experimental total pressure profile specified at the inlet section, the computed static 
pressure profiles are compared with the measurements in figure 14. The plots refer to the static pressure 
distribution in the section corners on the side wall and the symmetry plane. It is evident that the way the 
outer viscosity is computed may play a significant role in the correct prediction of the pressure distribution. 
Wall treatment 1, in which only the closest wall is considered, and 2, in which the outer turbulent viscosity 
is weighted according to the inverse of the wall distance, do not show large changes even if the two 
approaches are considerably different. For both techniques, the agreement with experimental results is 
fairly good on the pressure side of both the side wall and the symmetry plane. The suction side on both 
planes shows that the static pressure is overestimated. This is probably due to the presence of computed 
secondary flows much stronger than the experimental ones. The computed pressure drop induced by the 
presence of the bend and the strong flow acceleration is smoother than the measured one: This phenomena 
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is caused by high turbulent viscosities that induce a heavy momentum diffusion in the cross-flow direction 
followed by a static pressure redistribution and growth of the boundary layer thickness. The agreement 
with measurements improves with the turbulent viscosities averaging given by wall treatment 3. No big 
changes between the multiple wall treatments are found for the pressure side of both the side wall and the 
symmetry plane where it was always possible to have accurate results. Differences start to appear on the 
suction side where multiple wall treatment 3 shows the closest agreement with experiments. Still the 
pressure minimum located at 5=2, where velocities on the cross section start to develop, is not captured. 
This test suggests that the averaging technique based on the Van Driest damping expression for the mixing 
length can give reasonable predictions. All of the computations with the Baldwin-Lomax model were 
performed with multiple wall treatment 3. 

The results of the mesh dependence tests for the Baldwin-Lomax model are summarized in figure 15, 
where the computed static pressure distribution in the four corners of the cross sections using different grids 
are compared with measurements. The pressure distribution profiles show that there are practically no 
differences between the predictions obtained with grids 1 and 3. This proves that an increment of the grid 
points number in the main flow direction does not produce any gain in terms of accuracy: This is strictly 
connected to the boundary layer resolution that is not improved by using grid 3. The growth of secondary 
flows is influenced by the low momentum regions located close to the wall, the correct simulation of which 
is not ensured by the two grids. Conversely, the implementation of grid 2 clearly improves the accuracy of 
the results. The static pressure profile on the suction side of the side wall is in better agreement with 
experiments than the pressure distributions given by grids 1 and 3. It is worthwhile observing that the use 
of a more refined grid in the cross flow direction shows that the computed pressure minimum on the 
suction side is correctly located, even if its value is still overestimated, whereas this local minimum, located 
approximately at 5=2, is completely lost with the other two grids. Nevertheless, the static pressure 
distribution on the symmetry plane appears to be weakly affected by grid refinement. No further 
investigation was performed by varying the cross-flow points expansion ratio. 

The flow simulation with the q-u> model required more tests since it was necessary to investigate both 
the dependence on the mesh refinement and on the inlet turbulence level. The use of grids with expansion 
ratios equal to 1.2 did not ensure significant improvements in results since not enough points were located 
close to the wall. In order to have a reasonable definition of the turbulent kinetic energy peak at the wall, 
the grid expansion ratio was fixed at 1.3. The results of the first set of tests are summarized in figure 16 in 
which, adopting the coarse grid 4, the inlet turbulence level was changed to verify its influence on the static 
pressure distribution. Even when the inlet turbulence level is decreasing from 5.0% to 0.1%, the computed 
profiles progressively approach the measurements, the predictions are still far from experiments for the 
suction side of both the symmetry plane and the side wall. This indicates the presence of a large 
momentum diffusion that is possibly caused by insufficient mesh refinement or too high turbulence level, or 
both. Figure 17 shows the results obtained with grid 5 and lower turbulence levels. Figures 16 and 17 show 
that the refined grid with the same turbulence level brings some improvement of the agreement with 
experiments, proving that the coarser grid was largely inadequate for this flow configuration (see, for 
instance, the 0.1% level). The growth of the secondary flows is evident in figure 17 where the static pressure 
distribution computed with the inviscid approach is compared with the results of the two-equation model 
obtained with different turbulence levels. 

The direct comparison of the turbulent calculations with the inviscid computation, in which the same 
pressure distribution is found on both the side wall and the symmetry plane because of the absence of 
secondary motions, shows large differences on the suction side only starting from 5=2 where experiments 
deviate from the inviscid solution. The 1% turbulence level, not reported in figure 17, appeared to bring 
still a too high momentum diffusion so that this level was progressively decreased from 0.5% to 0.1%. The 
static pressure distribution on the pressure wall is mainly driven by convective phenomena, whereas the 
distribution on the suction wall is largely influenced by diffusion processes: this is the reason why the 
pressure side distribution is always reproduced well and is nearly independent of the inlet turbulence level. 
The final result obtained with grid 5 and a 0.1% turbulence level shows a fairly good agreement with 
experiments on both the suction and pressure sides and seems to reproduce quite correctly the location and 
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the influence of secondary flows. 

The predicted pressure distributions obtained by the Baldwin-Lomax model with wall treatment 3 
and the q-w model are very similar, but the two-equation model proved to be marginally more accurate 
especially on account of the static pressure distribution on the suction side. Figure 18 shows that the two 
models predict approximately the same static pressure pattern on the symmetry plane with nearly the same 
pressure drop due to the acceleration of the flow. Still, from figure 19, where the Mach number isolines are 
shown on the same symmetry plane, it is possible to observe that the zero equation model predicts a 
slightly thicker boundary layer than the one predicted with the q-w. These differences start immediately 
after the inlet section and become more evident as the exit section is approached. This indicates a different 
turbulent viscosity distribution in the wall region. The Baldwin-Lomax model was in fact found to predict 
a sharper growth of the turbulent viscosity in the wall region: still the two models gave comparable values 
of the turbulent viscosity in the flow core. The differences in the boundary layer thickness are evident in 
figure 20 where the velocities in the exit section of the channel are plotted. Although both the turbulence 
models show the same flow pattern, the two-equation model predicts the center of the secondary 
recirculation closer to the wall than that predicted by the zero-equation model in which the location of the 
center appeared to be in better agreement with experiments. This confirms a weaker cross-flow momentum 
diffusion given by the two-equation model as compared with the zero-equation model. This is caused by the 
aforementioned differences in the turbulent viscosities. It is remarkable that both formulations predict a 
small secondary vortex in the wall corner of the suction and pressure walls. 

An interesting qualitative comparison of the predicted and measured secondary flows is given by 
figure 21. Figure 21(a) shows the experimental flow visualization by injecting a smoke filament inside the 
boundary layer close to the side wall at the inlet section of a reduced model of the channel. The low- 
momentum particles located well within the boundary layer are bent towards the suction wall by the 
pressure gradients, but the high-momentum particles exhibit a weaker turning. Figure 21(b) shows the 
secondary flow predicted by the zero equation model, and figure 21(c) refers to the results obtained with 
the q-w model. From this comparison it is evident that the zero-equation model predicts a smoother 
turning of the velocities which is in slightly better agreement with the experimental picture than the q-w 
model. 

A typical total pressure loss distribution on the channel exit section is given in figure 22. The results 
computed with the Baldwin-Lomax model and grids 1 and 2 show the much thicker boundary layer 
obtained with the coarser grid. The refined grid ensures a correct description of secondary flows together 
with the location and magnitude of the total pressure losses. 

Concerning computational details, with the present research version of the code, 6000 sec were 
necessary to perform 2500 iterations with the q-w turbulence model with a 63.550-point grid on a Cray Y- 
MP supercomputer, which ensured an overall residual of the order of 10~ 6 , whereas a reduction of 
approximately 30% in CPU time could be obtained by using the zero-equation model. As mentioned by 
Coakley (ref. 15), the q-w model was found to be remarkably insensitive to the initial field specified to 
start the calculations for q and u>. In the present tests it has always been possible to start the model with a 
flat distribution of the turbulence quantities. For the viscous calculations it was found to be good practice 
to keep CFL<10 together with 1/32 as the weight for the fourth order artificial damping. The total 
absence of shocks allowed the second order artificial damping to be eliminated in all calculations. Moreover, 
the aforementioned wall treatment of the fourth order damping ensured inlet-outlet mass errors of the order 
of 0.5%. 


CONCLUDING REMARKS 

The scalar form of the approximate factorization coupled with a turbulence model proved to be 
suitable for solving turbulent internal flows where diffusive terms play a dominant role. The introduction of 
an approximate treatment of the diffusive terms proved to increase convergence of the algorithm for 
internal flow configurations. Nevertheless, the implicit treatment of these terms must be tested for a wider 
range of geometries and Reynolds numbers. The influence of Mach number on convergence rate needs to be 
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investigated especially for the proposed approximate implicit treatment of diffusive terms based on a space 
derivative of the fluid density: New tests are currently being performed. 

For the Stanitz elbow geometry the Baldwin-Lomax turbulence model proved to give results in 
acceptable agreement with experiments provided that the presence of multiple boundaries is properly 
accounted for. The low-Re q-w two-equation model version investigated here proved to be suitable for 
three-dimensional computations and gave a satisfactory description of the flow field with a manageable 
number of grid points and only a small increase of computational time with respect to the algebraic model. 
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Figure 7.b. S-duct (laminar flow regime ): Velocity profiles on midspan. 



Figure 7 .c. S-duct (laminar flow regime): Secondary velocity at bend exit. 






Figure 10. a. Stanitz elbow: Inviscid pressure 
isolines. 


Figure 10. b. Stanitz elbow: Inviscid Mach number 
isolines. 
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Figure 11. Stanitz elbow: Inviscid C p distribution. 
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Figure 12. Stanitz elbow: Convergence history. 
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Figure 14. Stanitz elbow: C p distribution with multiple wall treatments by Baldwin-Lomax model. 
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Figure 15. Stanitz elbow: C p dependence on grid points by Baldwin-Lomax model. 
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Figure 16. Stanitz elbow: C p dependence on turbulence level by q-ui model (coarse grid). 
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Figure 17. Stanitz elbow: C p dependence on turbulence level by q-u> model (refined grid). 
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Figure 18. Stanitz elbow: Static pressure isolines on symmetry plane. 
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Figure 19. Stanitz elbow: Mach number isolines on symmetry plane. 
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(a) Baldwin-Lomax (b) q -w 

Figure 20. Stanitz elbow: Velocity vectors at exit section. 
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Figure 22. Stanitz elbow: Total pressure losses. 
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